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ABSTRACT 

Contemporary pulsar timing experiments have reached a sensitivity level where 
systematic errors introduced by existing analysis procedures are limiting the achiev- 
able science. We have developed tempo2, a new pulsar timing package that contains 
propagation and other relevant effects implemented at the 1 ns level of precision (a 
factor of ^ 100 more precise than previously obtainable). In contrast with earlier tim- 
ing packages, tempo2 is compliant with the general relativistic framework of the lAU 
1991 and 2000 resolutions and hence uses the International Celestial Reference System, 
Barycentric Coordinate Time and up-to-date precession, nutation and polar motion 
models. Tempo2 provides a generic and extensible set of tools to aid in the analysis 
and visualisation of pulsar timing data. We provide an overview of the timing model, 
its accuracy and differences relative to earlier work. We also present a new scheme 
for predictive use of the timing model that removes existing processing artifacts by 
properly modelling the frequency dependence of pulse phase. 
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1 INTRODUCTION 

Pulsar timing observations have produced some of the most 
exciting results in pulsar astronomy and indeed in all of as- 
tronomy. For instance, such results have included the first 
detection of extra-solar planets (Wolszczan & Frail, 1992), 
stringent tests of the general theory of relativity (e.g. Stairs 
2003), revealed dispersion measure variations due to the in- 
terstellar medium (e.g. Backer et al. 1993), pulsar proper 
motions (e.g. Hobbs et al. 2004) and irregularities in the 
spin-down of pulsars (e.g. Lyne 1999). Pulsar timing is now 
being used to verify terrestrial time standards and the Solar 
System ephemeris and in searches for gravitational radia- 
tion (see for example, Foster & Backer 1990 and Jenet et al. 
2005). 

An overview of pulsar timing has been given by nu- 
merous authors (for example, Manchester & Taylor 1977, 
Backer & Hellings 1986, Lyne & Smith 1998 and Lorimer & 
Kramer 2005). In brief, the arrival times of pulses (TO As) 
are measured at a radio observatory for a particular pulsar 
over many years. These TO As need adjustment so that they 
represent arrival times in an inertial reference frame. This 
is accomplished by transforming each measured arrival time 
to an arrival time in the reference frame of the pulsar by 
first calculating arrival times at the Solar System barycen- 
tre (SSB) and then, if necessary, including additional terms 
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required to model the pulsar's orbital motion. A model of 
the pulsar's spin-down behaviour, the "timing model" or 
"timing ephemeris", is fitted to these arrival times. If sig- 
nificant systematic deviations are seen when calculating the 
differences between the actual arrival times and the best-fit 
model arrival times (known as timing residuals) then it is 
clear that the model is not fully describing the true pulsar 
parameters; a positive residual corresponds to the pulse ar- 
riving later than predicted. Such discrepancies can be due 
to many effects including unmodelled binary companions or 
binary parameters, irregularities in the spin-down of the pul- 
sar, or poor estimation of the astrometric or rotational pa- 
rameters. For instance, an incorrect estimate of the pulsar's 
position or its proper motion leads to a poor determination 
of the barycentric arrival times, which will produce a sinu- 
soidal feature in the timing residuals. This timing technique 
therefore allows pulsar parameters to be measured extremely 
precisely; the precision improves with longer data sets and 
more accurate TOA measurements. 

Both the conversion from the measured TOAs to 
barycentric arrival times and the model fitting required to 
obtain precise pulsar parameters are complex and can only 
be carried out within a computer program. Programs such 
as PSRTIME at Jodrell Bank observatory, timapr at Bonn, 
ANTIOPE at Nancay and CPHAS at Hartebeesthoek obser- 
vatories have already been developed. However, the most 
widely used and best known package is tempo which has 
been maintained and distributed by Princeton University 
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and the Australia Telescope National Facility^ . This package 
is extremely powerful, but the algorithms implemented are 
poorly documented and only provide a timing precision of 
~100ns. Recent high precision timing experiments produce 
root-mean-squarc (rms) residuals of this order and there- 
fore such results are systematically affected by inaccuracies 
in the TEMPO algorithms. A further limitation of TEMPO is 
that it can only be used to analyse one pulsar at a time. In 
order to study the recently discovered double-pulsar system 
(Lync ct al. 2004), to search for gravitational waves or to 
look for irregularities in terrestrial time standards, it is ad- 
vantageous to analyse multiple pulsars simultaneously. We 
have developed a now package, known as tempo2, which is 
based on the original tempo (hereafter called tempoI), but 
has a significant number of new and improved features. 

The aim of this paper is not to provide a user manual, 
but rather to 1) give a succinct description of the algorithms 
implemented, 2) highlight features that are not available in 
existing timing packages and 3) describe the accuracy of 
TEMP02. Full documentation and download instructions for 
TEMP02 can be obtained from our web site^. Details of the 
algorithms used in TEMP02 in order to achieve accuracies of 
1 ns will be presented in Paper II of this scries. Methods to 
simulate the effects of gravitational waves on pulsar timing 
data and utilities to place limits on the existence of a grav- 
itational wave background will be described in Paper III. 

In §2 we describe real and simulated pulse arrival times 
used for testing and demonstrating the various features of 
TEMP02. §3 provides a description of the conversion between 
site arrival times to arrival times in the pulsar frame through 
the use of clock correction files, propagation delays and a 
planetary ephemeris. The fitting algorithms implemented in 
TEMP02 for single datasets arc described in §4. §5 describes 
analysis methods for the fitted parameters and their uncer- 
tainties and §6 contains information on TEMP02 routines 
to study the resulting timing residuals. Tempo2 provides a 
predictive facility which is described in §7. 



2 REAL AND SIMULATED PULSE ARRIVAL 
TIMES 

The TEMP02 software is based around 1) an "engine" that 
calculates the barycentric arrival times, forms the timing 
residuals and carries out the weighted least-squares fit and 
2) "plug-ins" that add to the functionality of TEMP02 and 
allow the results to be analysed and presented in a user- 
friendly form. For instance, a plug-in is available to plot the 
timing residuals of multiple pulsars simultaneouslj', another 
to determine the power spectrum of the residuals and an- 
other to graph the clock corrections that tempo2 is applying 
to the measured arrival times. A full listing of the currently 
available plug-ins is provided in the Appendix. 

It is now common to combine TOAs obtained at dif- 
ferent observatories with different back-end systems and re- 
ceivers. These almost invariably give rise to a constant off- 
set or "jump" , between each set of TOAs. Tempo2 can fit 
for such jumps between observations at different telescopes, 

^ http://www.atnf.csiro.au/research/pulsar/tempo 
^ http://www.atnf.csiro.au/research/pulsar/tempo2 
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Figure 1. The timing residuals for PSRs J0437-4715, 
J1022-I-1001 and J1909-3744 observed as part of the Parkes Pul- 
sar Timing Array project. This figure was produced using the 
SPLK interface to tempo2. 



Table 1. A selection of the pulsars observed in the Parkes Timing 
Array Project. We list the pulsars' names, pulse frequencies (u), 
observing spans, numbers of TOAs (Nxoa)) observing frequencies 
(/) and the post-fit rms residuals (rms). 



PSR u Span Nxqa f rms 

(Hz) (d) (MHz) (us) 

J0437-4715 173.6879 761 6382 1340 0.49 

J1022+1001 60.7794 600 142 3100/1400/685 3.7 

J1909-3744 339.3157 542 56 3100 0.35 



with different observing frequencies or back-end systems, be- 
tween a range of dates, or on any other given parameter. This 
is made possible by a new free format for the measured pulse 
arrival times. This free format allows additional flags provid- 
ing user-definable parameters such as the backend system, 
the observation length or the observation bandwidth. Using 
the graphical plug-in features it is also possible to plot the 
pre- and post- fit timing residuals versus time or other pa- 
rameters such as the observation length, parallactic angle or 
attenuation settings. 

It is essential to test the algorithms implemented within 
TEMP02 with precise TOAs. We have selected three pul- 
sars, listed in Table 1 that have been observed for the 
Parkes Pulsar Timing Array project (PPTA)^ which aims 
to detect gravitational waves by looking for correlated sig- 
natures in the timing residuals of multiple pulsars (Jenet 
et al. 2005). Timing residuals for these pulsars are shown 
in Figure 1. PSRs J0437-4715 and J1909-3744 provide 
high quality timing observations at the sub-500ns level. 
PSR J1022-M001 has an ecliptic latitude of -0.06° and 
hence, the TOAs are affected by Solar System dispersion 
and Shapiro delays. Full details of the project and observ- 
ing details will be described in a later paper. 

For more detailed tests we use a tempo2 plug-in capable 

^ http: / /www.atnf.csiro.au/research/pulsar/ppta 
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Figure 2. Tempo2 simulation of "white" (upper panel) and "red" 
(lower panel) timing residuals. The lower panel contains both red 
timing noise simulated using a steep power-law spectrum and a 
small glitch event at MJD 51900. This figure was produced using 
the FAKE and Splk interfaces to tempo2. 

of simulating pulsar timing residuals in the presence of red 
noise or with glitch events (see Figure 2). These TO As are 
determined by repeatedly forming the pulsar timing resid- 
uals and then subtracting these residuals from the TOAs 
until the TOAs exactly match the timing model provided. 
The simulated residuals are then output after the addition 
of "white" (Gaussian) and/or "red" noise (modelled by sum- 
ming many sinusoids with random phase, but with ampli- 
tudes given by the requested power- law spectrum). 



3 FORMING THE PULSE EMISSION TIME 

The timing procedure starts by converting the measured 
topocentric TOAs to the pulse emission time in the pul- 
sar frame ignoring the frequency independent propagation 
delay from the pulsar to the SSB. Full details of this trans- 
formation, its accuracy and differences relative to tempoI 
will be described in Paper II. Here we summarise the trans- 
formation as 

At = Ac -I- Aa -I- Aeq -I- Arq + Asq - D/f + Ayp + Ab (1) 

where Ac contains various clock corrections (see §3.1), Aa 
the atmospheric propagation delays (§3.2), Aeq the Solar 
System Einstein delay (§3.3), Arq the Solar System Roemer 
delay (§3.4), Asg the Solar System Shapiro delay (§3.4.1), 
D/ p models the dispersive component of the light travel 
time (§3.5), Avp describes the excess vacuum propagation 
delay due to secular motion (§3.6) and Ab contains terms 
that describe any orbital motion (§5.3). In Tabic 2 wc list 
various effects that must be taken into account when forming 
barycentric arrival times from the observed TOAs. The table 
also provides a typical value or range for the magnitude of 
each effect and whether or not it is included in tempoI. 

3.1 Clock corrections 

The TOAs provided to tempo2 are recorded against local 
observatory clocks. Such clocks are typically derived from a 
precision frequency standard with good short-term stability. 



such as a hydrogen maser. On longer time scales (months to 
years) these clocks deviate significantly from uniformity and 
are therefore unsuitable for precision pulsar timing. How- 
ever, it is generally possible to remove these errors down to 
the precision provided by the best available terrestrial time 
scale through the apphcation of corrections derived from 
monitoring the offsets between pairs of clocks. For example, 
the PPTA pulsars are observed at the Parkes observatory 
where the offset between the observatory 1 pulse-per-second 
signal (derived from a hydrogen maser) is compared both 
to the clock signal broadcast by Global Positioning Sys- 
tem (GPS) satellites and by common-view GPS monitor- 
ing to the Australian national time scale UTC(AUS), main- 
tained by the National Measurement Institute. The Bureau 
International des Poids et Mesures (BIPM) in turn pub- 
lishes a monthly bulletin (Circular T) tabulating offsets be- 
tween various clock pairs. Using Circular T, measurements 
can be referred from an intermediate clock (e.g. UTC(AUS) 
or GPS time) to Universal Coordinated Time (UTC). UTC 
is a timescalc formed through the weighting of data from 
an ensemble of atomic clocks from around the world. This 
in turn is related to Temps Atomique International (TAI) 
by an integer number of "leap" seconds, which are inserted 
to maintain approximate synchrony between UTC and the 
irregular rotation of the Earth (these are announced in Bul- 
letin C of the International Earth Rotation Service). TAI is 
the most stable long-term time scale available in near real- 
time. 

The ultimate aim of the clock correction process is to 
transform measurements into the Geocentric Celestial Ref- 
erence System (GCRS), for which the coordinate time is 
denoted TCG, expressed in units of the SI second. Owing to 
their gravitational and rotational energy, terrestrial atomic 
clocks made to approximate the SI second do not run at the 
same rate as TCG. Instead, these clocks are used to define 
realizations of a timoscale known as Terrestrial Time (TT), 
which differs from TCG by a constant rate in such a way 
that its unit corresponds to the SI second on the surface of 
the geoid. One possible realization of TT is obtained directly 
from TAI: 

TT(TAI) = TAI + 32.184 s, (2) 

however TAI has instabilities and inaccuracies for which cor- 
rections frequently become available at a later date. The 
best available stability is currently provided by the retroac- 
tive timescales published by the BIPM (Guinot 1988, Petit 
2003), the most recent of which is denoted TT(BIPM04). 

The TEMP02 framework for handling clock corrections 
was designed with maxirrmm flexibility in mind, with the 
possibility of processing data sots with a heterogeneous col- 
lection of different observatories, clocks, and clock correction 
paths. The scheme is beised around a database of ASCII files 
tabulating the offsets between named pairs of clocks. Given 
the name of the clock against wfiicfi a TOA is measured, and 
the name of the realisation of TT to which it should be trans- 
formed, corrections can be applied based upon a manually or 
automatically determined sequence derived from linear in- 
terpolation of values from files found in the database. Step 
changes such as leap seconds are also possible. An ancillary 
suite of programs allows for the production of tempo2 for- 
mat files from external data sources such as Circular T, and 
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Table 2. Corrections and their typical sizes for phenomena included in tempo2. 



Correction 


Typical value/range 


TEMPOl 


Ob*iprvfltorv clnrW to ' 1 " 1 ' 




Y 




10 ns 


N 






N 


lAU precession/nutation 


~5ns 


N'^ 


Polar motion 


60 ns 


N 


AUTl 


IflS 


Y 


Einstein delay 


1.6 ms 


Y 


Roemer delay 


500 s 


Y 


Shapiro delay due to Sun 


112 //s 


Y 


Shapiro delay due to Venus 


0.5 ns 


N 


Shapiro delay due to Jupiter 


180 ns 


N 


Shapiro delay due to Saturn 


58 ns 


N 


Shapiro delay due to Uranus 


10 ns 


N 


Shapiro delay due to Neptune 


12 ns 


N 


Second order Solar Shapiro delay 


9 ns 


N 


Interplanetary medium dispersion delay 


100 ns' 


Y 


Interstellar medium dispersion delay 


~ls' 


Y 



" earlier precession/nutation model implemented 
observing frequency and pulsar dependent, typical value for 1400 MHz listed. 



provides capabilities for averaging, resampling and various 
analytic procedures for assessing the quality of data present. 

3.2 Atmospheric propagation delays 

The group velocity of radio waves in the atmosphere dif- 
fers from the vacuum speed of light. Rcfractivity is induced 
both by the ionised fraction of the atmosphere (mainly in 
the ionosphere) and the neutral fraction (mainly in the tro- 
posphere). The tropospheric propagation delay can be sep- 
arated into the so-called "hydrostatic" and "wet" compo- 
nents (see Paper II). For the highest timing precision, it is 
possible to provide tempo2 with a tabulated list of surface 
atmospheric pressure measured at an observatory for the cal- 
culation of the hydrostatic delay which will be of the order 
of 10 ns. If atmospheric pressure data arc unavailable then 
TEMP02 can, if required, use a canonical value of one stan- 
dard atmosphere. This assumption results in errors of the or- 
der of 1.5 ns. In Figure 3 we show computed hydrostatic tro- 
pospheric delays for simulated TOAs for PSR J1022+1001, 
assuming a constant surface atmospheric pressure and a ±5- 
h hour-angle range. Diurnal variations arise due to the de- 
pendence of atmospheric path length on source elevation (in 
the simulated observations the elevation varies from 6 to 46 
degrees) . 

The wet component of the tropospheric propagation de- 
lay (the zenith wet delay, ZWD) is highly variable and can- 
not be predicted accurately. If no tabulated ZWD informa- 
tion is available the effect is neglected, otherwise tabulated 
data may be used. With a typical excess zenith path length 
of 100-400 mm, error is incurred at the level of approxi- 
mately 1.5 ns. 



3.3 Einstein delay 

The Einstein delay (Damour & Deruelle 1986) quantifies 
the change in arrival times due to variations in clocks at 
the observatory and the SSB due to changes in the gravitar 
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Figure 3. The computed hydrostatic tropospheric delay for sim- 
ulated pulse times of arrival, assuming a constant surface atmo- 
spheric pressure. 



tional potential of the Earth and the Earth's motion. lAU 
resolution A4 (1991) recommends the use of barycentric co- 
ordinate time (TCB) which differs from TT both in mean 
rate and in periodic and quasi-periodic terms. By default, 
this is the coordinate time in which arrival times are speci- 
fied in TEMP02. Prior to the definition of TCB, the recom- 
mended barycentric coordinate time was Barycentric Dy- 
namical time (TDB) which was implemented in TEMPOl. In 
addition to being physically unrealisable (Standish 1998), 
TDB values are not physical coordinate times, but rather 
values of a variable related to time by a dimensionless scale 
factor (Klioner 2005). If these values are taken as barycen- 
tric coordinate times of arrival, as has been common practice 
in the past, then the scaling factor is effectively transferred 
from the value to the units. Therefore, although site arrival 
times are referred to TT, which is defined in terms of the 
SI second, TDB barycentric arrival "time" intervals, and 
in fact, the numerical values of all parameters inferred with 
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pulsar timing on the basis of TDB TOAs are effectively mea- 
sured in units that differ subtly from their SI counterparts. 

As a result, all catalogued parameters measured using 
TDB TOAs (e.g. those from tempoI) must be multiplied 
by 

K = l + (1.55051979154 x 10"® ± 3 x 10"") (3) 

(Irwin & Fukushima 1999). This can be a large effect and 
timing models created using tempo 1 need to be modified 
before being used by tempo2. The n'th frequency derivative 
scales as if and the orbital period and semi-major axis 

all scale as K. The epochs of periastron, period, position 
and dispersion measure all scale as K in their offset from 
the common epoch of Modified Julian Day 43144.0003725 
(Irwin & Fukushima 1999). For instance, the modification in 
pulse frequency produces a slope of 0.5 syr"^ in the timing 
residuals which, for millisecond pulsars, will lead to phase 
coherence being lost over even short data spans. Using the 
TRANSFORM plug-in, TEMP02 provides an interface that can 
be used to convert old parameters into the new system. We 
also emphasise that, because of the significant differences 
between the TDB and TCB, for all published timing models 
the coordinate frame used must be clearly specified. 

3.4 Roemer delay 

The Roemer delay is the vacuum light travel time between 
the pulse arriving at the observatory and the equivalent ar- 
rival time at the SSB. In tempo2 this is calculated by de- 
termining the time delay between a pulse arriving at the 
observatory and at the Earth's centre and, with the aid of 
a Solar System ephemeris, from the Earth's centre to the 
SSB. 

The coordinates of the pulsar are known, either from 
telescope pointing, interferometry or pulsar timing, and are 
normally measured in the International Celestial Reference 
System (ICRS)''. The required transformation between the 
ICRS and the International Terrestrial Reference Frame 
(ITRF), within which observatory positions are determined, 
depends on precession, nutation, polar motion and Earth 
rotation. The worst-case timing offset resulting from a posi- 
tional error of is given by A^J?0 /c, i.e. 



1 ns 9.7 mas ' ^ ' 

Through its omission of polar motion amounting to up 
to ± 300 mas (corresponding to ±30 ns) and also through 
the use of lAU 1976 precession (Lieske et al. 1977) and lAU 
1980 nutation (Seiber 1982) models which are in error at the 
50 mas (5 ns) level, the tempo 1 software introduces errors in 
the timing model that are significant at contemporary levels 
of timing precision. In tempo2, polar motion is corrected 
using the values published in the C04 series of Earth Orien- 
tation Parameters (EOF) of the International Earth Rota- 
tion Service (lERS). The lERS also provides the difference 

* In the case of positions obtained by pulsar timing, this is only 
true if the reference frame of the Solar System ephemeris is tied 
to the ICRS, e.g. by using the DE405 planetary ephemeris. The 
DE200 ephemeris is offset from the ICRS by ~14 mas (Folkner 
et al. 1994) , yielding a potentially significant error in the transfer 
to the geocentre if such a position is used; see Paper II. 



between the observed precession and nutation and that pre- 
dicted by the lAU 1976 and 1980 models. However, follow- 
ing the recommendations of lAU Resolutions adopted at the 
24th General Assembly, we adopt the lAU 2000 precession- 
nutation model which provides sufficiently accurate predic- 
tions. Specifically, tempo2 uses the truncated 2000B model 
(McCarthy & Luzum 2003) which is accurate to 1 mas 
(0.1ns). Figure 4 shows the differences between the Solar 
System Roemer delay computed using TEMPOl and TEMP02 
using simulated observations of PSR ,11022+1001. Differ- 
ences in the model (mainly due to polar motion) introduce 
an error in the assumed observatory position, which appears 
as a diurnal timing term which is modulated by the yearly 
and 435-d periodicities of the polar motion. 

The third component in the transformation of the pul- 
sar position to the ITRS is the Earth rotation angle which 
is a linear function of the time scale known as UTl. This 
is computed by tempo2 using the offset between UTC and 
UTl as provided in the C04 EOP series. 

The choice of Solar System ephemeris for determining 
the position of the SSB with respect to the Earth can have 
significant effects on the calculated timing residuals. Until 
recently the JPL DE200 model (Standish 1990), which is 
based upon the dynamical equator and equinox of J2000, 
was the most widely used. More recently, the JPL DE405 
model has been developed^ which, in contrast to the DE200 
model, is aligned with the ICRS (Standish 1998). The DE405 
model includes the planets, the Earth's Moon and 300 as- 
teroids. In the left panel of Figure 5 we plot the difference 
between residuals obtained using the DE405 and DE200 
models after the subtraction of an annual sinusoid (cor- 
responding to a position error) and quadratic term (cor- 
responding to the spin-frequency and its first derivative) 
for simulated observations of PSR J1909-3744. The right 
panel contains the timing residuals after fitting for five fre- 
quency derivative terms, the orbital period, epoch of peri- 
astron, proper motion and parallax. The use of the DE200 
model leads to an incorrect measurement of the proper mo- 
tion in right ascension by —0.2775(3) masyr"'^, in declina- 
tion by —0.037(1) mas yr"^ and parallax by —0.045(5) mas 
over this simulated, regularly sampled data span of 14 years. 
Splaver et al. (2005) also reported significant deviations be- 
tween residuals for PSR J1022-I-1001 using these Solar Sys- 
tem ephemerides during the years 1998-1999 which they ex- 
plained by the new ephemeris incorporating improved mea- 
surements of the outer planet masses. Although tempo2 can 
access any of the JPL planetary ephemerides, we currently 
recommend that the DE405 model be used for any high pre- 
cision analysis of pulse arrival times. 



3.4.1 Shapiro delay 

To make an accurate determination of the arrival time at 
the barycentre it is also necessary to include the Shapiro 
delay due to Solar System objects (most notably the Sun) 
which accounts for the time delay caused by the passage of 
the pulse through large gravitational fields (Shapiro 1964). 
Table 2 shows the maximum variation in Shapiro delay for a 



® ftp://ssd.jpl.nasa.gov/pub/eph/export/DE405/de405.iom/ 
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Figure 4. Differences in tlie Solar System Roemer delay computed using current lAU precession-nutation models and including polar 
motion, versus the algorithm of TBMPOl for simulated TOAs from PSR J1022-I-1001. The diurnal timing term is shown in the leftmost 
plot which is modulated by the yearly and 435-d periodicities of the polar motion (right) which, in turn, beat with a ~6-year period. 
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Figure 5. Comparison between timing residuals obtained using the DE200 Solar System ephemeris and the DE40-5 cplicmeris. In the 
left plot, terms corresponding to a pulsar position error, spin-frequency and its first derivative have been subtracted. In the right plot, 
terms corresponding to the above and higher frequency derivatives, orbital period, epoch of periastron, proper motion and parallax have 
also been removed. The vertical lines are spaced at lyear intervals. These plots were created using the plk plug-in for tempo2. 



selection of Solar System bodies. Tempo2 includes all bod- 
ies for which the maximum variation is greater than 0.1 ns. 
Figure 6 shows the variations in the Shapiro delay due to 
Jupiter for PSR J1022-M001®. The effect of the Shapiro de- 
lay due to the Sun can be clearly seen in the observations of 
PSR J1022-hl001 which has an ecliptic latitude of -0.06°. In 
Figure 7a we plot the pulse timing residuals after fitting for 
the pulsar's parameters and the Solar System Shapiro delay. 
Figure 7b shows the resulting timing residuals if the best-fit 
parameters are used, but the Solar System Shapiro delay is 
not calculated when forming the barycentric arrival times. 
On MJD 53143, this pulsar passed within 5° of Jupiter. How- 
ever, the additional Shapiro delay due to Jupiter is not de- 
tectable with our current data. 
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Figure 6. The additional time delay from the Shapiro delay due 
to Jupiter for PSR J1022-I-1001. 



® The Shapiro delay as characterised by Damour & DereuUe 
(1986) can be negative. However, as the zero-order time of ar- 
rival of the pulses is arbitrary, a constant offset can be added to 
the Shapiro delay calculation. 



3.5 Frequency-dependent peirameters 

Tempo2 provides the ability to fit for delays which are de- 
pendent upon the observing frequency; see Table 3. For in- 
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Figure 7. The timing residuals, in /is, for PSR J1022+1001, a) after fitting for the pulsar's parameters and b) without removal of the 
Solar System Shapiro delay. This plot was created using the plk plug-in for tempo2 (note, the original plk plotting package incorrectly 
plotted the uncertainties on the residuals; the errors were a factor of two too small). 



Table 3. Tempo2 parameters relevant to frequency-dependent offsets. 



Parameter Description Symbol 

DM, DM1 . . . The dispersion measure and its derivatives DM, DM . . . 

DMEPOCH The epoch of the disperison measure (MJD) to 

FDDI Index for frequency dependent delay C, 

FDDC Scale for frequency dependent delay kf 



stance, dispersion measure (DM) delays are oc whereas 
delays caused by refractive and diffractive effects are oc 
(e.g. Foster k. Cordes 1990). Tempo2 allows fitting for a pa- 
rameter that is oc where C, is defined by the user and is 
not restricted to integral values. Wc emphasise that in or- 
der to obtain absolute values for these frequency-dependent 
terms it is necessary to obtain TOAs using aligned stan- 
dard templates. In practice, true absolute alignment is im- 
possible because of profile shape evolution with frequency, 
so frequency-dependent parameters are always relative at 
some level. 

Although DM values are commonly published, the di- 
rectly measurable parameter from pulsar timing observa- 
tions is D, the dispersion constant, where 

DM = D/ku. (5) 

If the effect of ions and magnetic fields in the interstellar 
medium are ignored, then 

A;d = (6) 

However, ions and magnetic fields introduce a rather uncer- 
tain correction of order a part in 10^ (Spitzer 1962), compa- 
rable to the uncertainty in some measured DM values (e.g. 
Phillips & Wolszczan 1992). Consequently, both tempoI 
and TEMP02 adopt a value of ko = 2.410 x 10"^'^ cm~^pcs 
(Manchester & Taylor 1977). It is also possible, in tempo2, 
to set fcD = 1 in order to measure the dispersion constant. 

Another dispersive component occurs in the Solar Sys- 
tem. The interplanetary medium is dominated by the So- 
lar wind and is approximated in TEMP02 with the electron 



density decreasing as an inverse square law from the cen- 
tre of the Sun (full details are provided in Paper II) with 
no being the electron density at the Earth. TempoI uses 
no = 9.961 cm~^. However, by default, TEMP02 uses a value 
of no = 4cm^'' which is more consistent with recent mea- 
surements (Issautier et al. 1998). Figure 8 shows the extra 
time delay added by tempo2 for simulated observations of 
PSR J1022+1001. As discussed in Paper II, this estimation 
of the extra time delay is poor as the true electron den- 
sity can vary dramatically. We therefore recommend that, 
for high-precision timing, tempo2 be provided with multi- 
ple frequency observations which allow the determination of 
the actual DM for each observation. 

As an example, the PPTA uses a dual-band receiver 
at 10 and 50 cm. The DM at any instant can therefore be 
measured to a precision of up to 1 x 10~* cm~^pc if the 
difference between observations at the two frequencies can 
be measured to 1 /is. Tempo2 can be run in a mode where 
simultaneous (or contemporaneous) observations at multiple 
frequencies are used to determine the current DM and to use 
that value in subsequent calculations. For times when there 
are no multi-frequency observations available, the DM can 
be estimated from a polynomial fit to measured DM values 
before and after the observation. 

3.6 Shklovskii effect and radial motion 

Pulsar timing measurements are affected by the secular mo- 
tion of the pulsar relative the SSB. In the past the secular 
terms involving this motion have been omitted from timing 
models, because they can be absorbed in alterations of other 
parameters. The four largest effects are the radial veloc- 
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Figure 8. The extra time-delay automatically added by tempo2 
to model the interplanetary medium for PSR J1022+1001 at an 
observing frequency of 1400 MHz. This figure was produced using 
the PAKE and delays plugins for tempo2. 
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Figure 9. Comparison between TEMPOl and TEMP02 timing 
residuals for the PSR J0437-4715 observations when tempo2 
is emulating tempoI. No clock corrections were applied to the 
TOAs. 



ity (affecting most spin and binary parameters; Damour & 
Deruelle 1986), the Shklovskii effect and radial acceleration 
(affecting the spin and orbital period derivatives; Shklovskii 
1970, Damour & Taylor 1991) and the mixing of radial ve- 
locity into the Shklovskii term (affecting the spin period sec- 
ond derivative; van Straten 2003). In contrast to tempo 1, 
TEMP02 takes the approach that these terms can be included 
in the timing model as long as steps are taken to ensure 
the model is sufficiently constrained. In this way, one may 
take into account what is known about the secular motion 
and distance (via, for example, its appearance in the annual 
proper motion and parallax terms) to provide correct mea- 
surements of the spin and orbital parameters, rather than 
measuring incorrect values and attempting to correct them 
post facto (e.g. Damour & Taylor 1991). Conversely, if one 
may safely assume that one of the affected spin or orbital 
parameters is zero, it may be held fixed at this value in order 
to obtain a direct measurement of the distance or velocity, 
rather than measuring incorrect spin and orbital parame- 
ter values and using these to infer the motion and distance 
indirectly (e.g. Bell & Bailes 1996, van Straten 2003). 



4 FITTING ROUTINES 

Tempo2 uses the derived time of emission and a given tim- 
ing model to form the i'th pre-fit timing residual: 



(j>i describes the time evolution of the pulse phase based on 
the model pulse frequency (u) and its derivatives iu addition 
to any glitch parameters. Ni is the nearest integer to 4>i- 
Paper II contains details for calculating <pi. 

Terms corresponding to offsets in model parameters are 
fitted to these residuals in order to improve the measure- 
ment of these parameters. By default, the entire procedure 
is repeated using the post-fit timing model in order to pro- 
duce accurate post-fit barycentric arrival times and residu- 
als. This is in contrast to tempo 1 which only obtains the 
barycentric arrival times once and predicts the expected 
post-fit timing residuals. This entire process often needs to 
be iterated until convergence is reached as the offsets made 



to model parameters are based on a linearised approxima- 
tion to the effects on the timing model (e.g. Damour & Deru- 
elle 1986). 

The fitting routines in TEMP02 are based on a linear 
singular-value decomposition, weighted least-squares algo- 
rithm (for example. Press et al. 1992)^ where 

is minimised. A*' is the number of observations and ai — 1 
for unweighted fits or set to the TOA uncertainty for the 
i'th observation for a standard x^-minimisation. If specified 
as Modified Julian dates, measured TOAs need to be ac- 
curate to better than 19 significant figures for 1 ns timing 
precision and spin frequencies are now routinely measured 
to 16 significant figures. In tempo2 all parameters are stored 
and all calculations are carried out with "long double" pre- 
cision that typically provides 12 bytes of storage (allowing 
18 significant digits) on PC-based systems and 16 bytes (33 
significant digits) on most other systems. In Figure 9 we 
plot the difference between the pre-fit timing residuals for 
PSR J0437— 4715 obtained using tempoI and tempo2. The 
observed trend that covers ~2ns over 2.1 years of observ- 
ing is due to TEMPO 1 not storing the pulse frequency with 
enough significant figures and will lead to tempoI introduc- 
ing systematic effects in the timing residuals over long time 
spans. 

It is also useful to compare the formal uncertainties 
described above with those obtained using a "bootstrap- 
ping" method (see, for instance. Wall & Jenkins 2003) which 
can produce more realistic parameter values and uncertain- 
ties when significant correlations between parameters are 
present. The bootstrapping method implemented in tempo2 
estimates the uncertainty on a parameter by 1) randomly se- 
lecting observations to produce a new dataset of the same 
length as the original (the observations are selected with 
replacement; i.e. in the new dataset some of the original ob- 

^ It is possible with plug-in capabilities to use a non-linear fitting 
algorithm, but this is not necessary for the routines described in 
this paper; more details will be provided in Paper III. 
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Figure 10. Histogram of the declination parameter obtained us- 
ing the bootstrap technique for PSR J1909-3744. The dashed 
lines indicate the Icr uncertainties obtained using this method 
and the dotted lines give the Icr uncertainties measured using the 
standard least-squares-fitting routine. 



servations will be omitted while others will be replicated) , 2) 
recalculating the parameter and 3) repeating as many times 
as possible*. The distribution of these parameters provides 
an estimate of the uncertainty (obtained from the standard 
deviation of the distribution) on the value of the parame- 
ter (taken as the mean of the distribution). In Table 4, we 
compare the values and uncertainties on the fitted param- 
eters using the formal least-squares fitting and the boot- 
strap technique with 1024 iterations for PSR J1909-3744. 
For this pulsar, with residuals that are not dominated by 
timing noise, the measured uncertainties are typically ~ 1.2 
times larger with the bootstrap technique than with the 
least-squares method. A histogram of the fitted declination 
is shown in Figure 10. 

Various interfaces exist that allow the user to study the 
actual fits being applied to the pre-fit timing residuals. For 
instance, in Figure 11a we plot the components of the fit 
when improving a pulsar's proper motion and parallax. In 
Figure lib we demonstrate the effects of poorly estimated 
binary parameters. The tempo2 software also provides the 
ability to fit one or more of the timing model parameters 
over short adjacent subsets of the data. For example, this al- 
lows the user to analyse dispersion measure variations (Fig- 
ure 12), search for glitch events and to confirm proper mo- 
tions. 

It is often difiicult, with sparse observations, to obtain 
an accurate timing model. This is often solved by making 
further observations of the pulsar, but with the gorilla 
plug-in, TEMP02 provides an alternative method. For exam- 
ple, 35 observations of PSR J0857— 4424 spread over seven 
years were obtained at the Parkes observatory, but a solu- 
tion producing phase-connected timing residuals could not 
be obtained. Gorilla provides a brute-force fitting tech- 
nique that obtains the pre-fit timing residuals over many 
millions of combinations of spin-frequency and its derivative 
within specified ranges. This method found the correct so- 
lution over a 1-year section of the PSR J0857— 4424 data 

* Note: the bootstrapping method implemented in tbmpoI does 
not re-fit for the parameters after randomly selecting observations 
and is therefore not a true bootstrapping technique. 
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Figure 1 1 . The components of a fit to update a pulsar's a) proper 
motion and parallax and b) longitude of periastron and projected 
semi-major axis of orbit (assuming that these are the only pa- 
rameters being updated). The line through the measured residu- 
als indicates the sum of the model components. These plots were 
created using the plk plug-in for tempo2. 

and this fit was extrapolated to phase-connect the entire 
seven years of observations without the necessity for further 
observations. 



5 THE POST-FIT PARAMETERS AND THEIR 
ERRORS 

5.1 Pulsar spin parameters 

Tempo2 provides the ability to fit for the pulsar's pulse fre- 
quency and an arbitrary number of spin-frequency deriva- 
tives. The full set of possible spin parameters available are 
listed in Table 5. Except for the very youngest pulsars, the 
frequency second and higher derivatives are not believed to 
represent the secular spin-down of the pulsar, but rather 
timing irregularities known as timing noise (e.g. Hobbs et 
al. 2004). For pulsars whose timing residuals are dominated 
by timing noise it is not possible to determine accurate 
positions, proper motions or dispersion measures without 
"whitening" the data while fitting the timing model. Tradi- 
tionally, this whitening procedure has been carried out by 
fitting multiple spin-frequency derivatives until the result- 
ing post-fit residuals are free of systematic structure. As 
described by Hobbs et al. (2004) this whitening technique is 
limited because 1) only low- frequency timing noise can be 
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Table 4. Comparison between standard least-squares (LS) parameters and uncertainties with those obtained using a bootstrapping (BS) 
technique for PSR J1909-3744. 



Parameter 


Value (Vls) 


Error (Els) 


Value (Vbs) 


Error (Ebs) 


Ebs/E^bs 


Right ascension (rad) 


5.016908214879 


3.5x10-11 


5.016908214880 


4.3x10-11 


1.2 


Declination (rad) 


-0.65863987098 


1.4x10-1° 


-0.65863987100 


2.1x10-1° 


1.5 


Pulse frequency (Hz) 


339.31568762926 


1.9x10-12 


339.31568762926 


2.4x10-12 


1.3 


Frequency derivative (s~^) 


-1.614873x10-15 


2.4x10-2° 


-1.614878x10-15 


2.9x10-2° 


1.2 


Orbital period (d) 


1.533449474188 


1.1x10-11 


1.533449474191 


1.4x10-11 


1.3 


Projected semi-major axis (It-s) 


1.897991295 


9.9x10-° 


1.897991295 


1.2x10-* 


1.2 


Epoch of periastron (MJD) 


52053.452 


0.021 


52053.443 


0.020 


0.95 


Eccentricity 


1.186x10"'^ 


9.5x10-° 


l.mxlO-'' 


9.6x10-° 


1.0 
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Figure 12. The dispersion measure variation for PSR B0458-I-46 

obtained by fitting for dispersion measure in 1 year sections of 
Jodrell Bank observatory data. This plot was produced using the 

STRIDEFIT plug-in for TEMP02. 



modelled without affecting the higher-frequency signatures 
of position errors and proper motions and 2) such whiten- 
ing is limited to polynomials of order 12 to prevent floating- 
point overflows. Hobbs et al. (2004) described a new method, 
based on the fitting of harmonically related sinusoids, that 
produces superior results in many cases. This siuusoidal fit- 
ting technique has been implemented into tempo2. In brief, 
TEMP02 obtains from the parameter file a fundamental fre- 
quency (uj) and the amplitudes of uh harmonically related 
sinusoids. If the fundamental frequency is not provided then 
it is derived from the observation span (Tspan) as 



Tempo2 subsequently subtracts 

Aii = ^ aj, 8m.{ku)M) + bk cos{kujAt) (10) 

k = l 

from the timing residuals where At represents the difference 
between an arrival time in the pulsar reference frame and 
the specified period epoch. Tempo2 can also fit for the co- 
efficients tti and bi and output these parameters as part of 
a new timing model. In contrast to the technique described 
by Hobbs et al. (2004), the default method implemented by 
TEMP02 simultaneously fits for the pulsar paxameters and 
for the sinusoids, i.e. this is not a pre-whitening technique. 



An example of such "whitening" of the data is shown 
in Figure 13. The timing residuals of PSR B1842-H4 ob- 
tained from the Jodrell Bank data archive are typical of 
those seen for non-recycled pulsars over many years of ob- 
servation. They are dominated by quasi-periodic structures 
that are well modelled using the sinusoidal modelling, but 
not by high-order polynomial terms. 

The pulse frequencies of many young pulsars have been 
observed to increase suddenly during a glitch event. An indi- 
vidual glitch can be characterised by the epoch of the glitch 
[tg), the phase increment at the glitch {A(f>g), a permanent 
pulse frequency increment (Avg), a permanent frequency 
derivative increment (Ai^g), a decaying pulse frequency in- 
crement {Avd) and the decay time constant (xd). The in- 
crease in pulse phase due to a glitch event at time tg is 
given by 

4>g = A<j>g + AtAvg + ^^AvgiAtf + TiAva{l - e-'^*/"<*)(ll) 

where At = t — tg is the time since the glitch event. Tempo2 
allows the user to fit for an axbitraxy number of glitch events 
in a single dataset and provides various plug-ins to aid in 
their detection. 

5.2 Astrometric parameters 

Tempo2 performs computations in the reference frame of 
the Solar System ephemeris which is approximately equa- 
torial. Since the location of the observatory is specified in 
the ITRS, and transformed to the ICRS according to the 
specifications of the lAU 2000 resolutions, for the best pos- 
sible accuracy a Solar System ephemeris that is aligned with 
the ICRS should be used. The JPL DE405 ephemeris is the 
most recent of the publically available JPL series and meets 
this criterion. The tempo2 equatorial astrometric param- 
eters ("RA", "DEC", etc.) strictly refer to the ephemeris 
frame, which, even in the case of DE405, is only tied to the 
ICRS within a finite uncertainty. The ICRS itself is mea- 
surably offset from both the Fundamental Katalog 5 (FK5; 
Fcissel & Mignard 1998) and the dynamical equator and 
equinox of J2000.0 (Hilton & Hohenkerk 2004). Seidelmann 
& Kovalevsky (2002) state that ICRS coordinates should be 
denoted "right ascension" and "declination" with no further 
qualification. Since the coordinates measured using pulsar 
timing are relative to the coordinate frame of the chosen 
planetary ephemeris, the latter needs to be specified when 
quoting fitted positions. Owing to the known offset between 
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Table 5. The spin, glitch and whitening parameters included in tempo2. 



Parameter 


Description 


Symbol 


FO, Fl . . . 


The pulse-frequency and its derivatives 




PEPOCH 


The epoch of the pulsc-frcqucncy measurement 


tp 


GLEP.k 


Glitch epoch (MJD) 


tg 


GLPH.k 


Glitch phase increment 


A4>g 


GLFO.k 


Glitch permanent pulse frequency increment (Hz) 


AVg 


GLFl.k 


Glitch permanent frequency derivative increment (s~^) 


AVg 


GLFODJc 


Glitch decaying pulse frequency increment (Hz) 




GLTDJc 


Glitch decay time constant (d) 


Td 


WAVE.OM 


Fundamental frequency of sinusoids for whitening (Hz) 




WAVEJc 


Amplitude of the sine and cosine terms for the k'th sinusoids 





1842+14 (rms = 43701.197 us) pre-fit 1842+14 (rms = 2311.638 us) post-fit 




-2000 2000 

MJD-49362.4 



Figure 13. a) The pre-fit timing residuals for PSR B1842-(-14 obtained from the Jodrell data archive, b) Whitening the timing residuals 
using 11 polynomial coefficients and c) whitening the timing residuals using sinusoids. 



J2000.0 and the ICRS we recommend that the common prac- 
tise of labelling pulsar timing coordinates as "J2000" be dis- 
continued. 

Tempo2 also accepts astromotric parameters specified 
in the ecliptic frame. For such parameters, tempo2 trans- 
forms all vectors into the ecliptic coordinate system by rotat- 
ing about the "x-eoiis" by the mean obliquity of the ecliptic, 
e, at the epoch J2000.0. By default, tempo2 uses the current 
best estimate of the mean obliquity of the ecliptic, eDE405 = 
84381.40578 arcsec (Harada & Fukushima 2004). This value 
is derived from a harmonic decomposition of the DE405 So- 
lar System ephemeris, and applies to the epoch J2000.0. This 



differs from the earlier estimate of 84381.412 arcsec used 
by TEMPO 1. Since timing pulsar positions are largely con- 
strained by the Solar System Roemer delay, error ellipses 
tend to be aligned with ecliptic latitude or longitude, mak- 
ing this basis a convenient choice when one coordinate is 
constrained more strongly than the other. For instance, fit- 
ting for the proper motion in equatorial coordinates with 
the PSR J1022-I-1001 dataset discussed in §2 gives a proper 
motion in right ascension fia = —148(55) mas yr~^ and in 
declination of ns = —335(143) mas yr~^ and the astromet- 
ric parameters are highly covariant in the fit. Fitting using 
ecliptic coordinates provides a precise measurement of the 
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proper motion in ecliptic longitude jj,x = —16.1(2) masyr"^ 
which can be used in 1-dimensional studies of pulsar veloc- 
ities (e.g. Hobbs et al. 2005) and a much poorer determina- 
tion in ecliptic latitude ng = —307(152) masyr"^. A full list 
of the astromctric parameters that can be used by TEMP02 
is listed in Table 6. 

5.3 Binary parameters 

For pulsars in binary systems, TEMP02 includes terms that 
describe the pulsar's orbital motion. Various timing mod- 
els are available. A full mathematical description of these 
models and their implementation in tempo2 is provided in 
Paper II. Here, we provide a summary. 

The main binary model implemented in tempo2 is re- 
ferred to as the 'T2' model and is based on the Damour & 
Deruelle (1986) model ('DD') implemented in TEMPOl. In 
contrast to the 'DD' model which is designed for a pulsar 
and a single companion, the new 'T2' model allows mul- 
tiple binary companions. Various other models were avail- 
able to TEMPOl and can be emulated using the 'T2' model. 
For instance, the 'DD' model is more general than the ear- 
lier Blandford & Teukolsky (1976) model ('BT'). By mak- 
ing various simplifying assumptions (see Paper II), the 'BT' 
model can be obtained from the 'T2' model. Recently it 
has been shown (Kramer & Wex, private communication) 
that the 'DD' model provides poor uncertainties for mea- 
surements of the orbital Shapiro delay when the orbital 
inclination i ~ 90°. In the 'T2' model (and the 'DDS' 
model) a new parameter x = — ln(l — sini) can be used 
for such edge-on binary systems. As shown by Kramer (pri- 
vate communication) tempo2 provides more reliable esti- 
mates of the uncertainty on the new parameter x (known 
as SHAPMAX as it relates to the maximum Shapiro de- 
lay for a near-circular orbit) than on sini for near edge- 
on binary systems. For instance, with our observations of 
PSR J1909— 3744, the 'DD' model gives a companion mass 
of rUc = 0.2061(16) M0 and sinz = 0.99820(8). The 'T2' 
model gives rric = 0.2063(17) Mq and x = 6.28(6) implying 
that sini = 0.99813(11). Figure 14 gives a plot indicat- 
ing the fitted parameters obtained using the 'T2' and 'DD' 
models. 

The 'T2' model, by default, returns theory-independent 
results. However, it is possible to assume that general rel- 
ativity applies, emulating the 'DDGR' model developed by 
Taylor (1987) and Taylor & Weisberg (1989) for tempo 1. 
This allows the determination of the mass of the pulsar and 
its companion. 

For wide-orbit binary systems, Kopcikin (1995) showed 
that orbital timing parallaxes are measurable. Kopeikin 
(1996) also described secular variations of orbital paxam- 
eters due to the system's proper motion. These terms are 
included as part of the 'T2' model, but require an estimate 
of the longitude of the ascending node of the binary's orbit, 

n. 

The PSR B1259— 63 system of a pulsar orbiting a 
Be star has proved difficult to model. Wang, Johnston & 
Manchester (2004) therefore developed the modified BT 
model that allowed for jumps in the Kcplcrian parameters at 
specified times (specifically, for this system, at periastron); 
such jumps are available in the 'T2' model. A more physical 
approach to the effect of non-point-mass companions was 




0.9978 0.998 0.9982 0.9984 



sin(inclination) 

Figure 14. contours in m2-sini parameter space for 
PSR J1909— 3744. The circle indicates the most-likely values ob- 
tained using the 'DD' model and the cross symbol for the 'T2' 
model. This plot was produced using the m2sini plug-in. 

developed by Wex (1998) through alterations in the secular 
behaviour of the longitude of periastron and the projected 
semi-major axis. These changes are also available as part of 
the 'T2' binary model. 

For orbits with small eccentricities, the 'T2' model, by 
default, produces highly covariant values for the epoch and 
the longitude of periastron. The 'ELLl' model of Wex (1998 
unpublished work; see Lange et al. 2001) was developed for 
pulsars in such orbits. Tempo2 provides a tool to convert 
between a timing model based on the BT model and one 
based on the ELLl model which is defined by the parameters 

EPSl = esinw (12) 

EPS2 = ecosw (13) 

TASC = To-^Pb (14) 

ZTT 

where e represents the orbital eccentricity, a; the longitude 
of periastron. To the epoch of periastron and Ph the orbital 
period. If these parameters are included in a 'T2' parameter 
file then the 'T2' model will emulate the 'ELLl' model. 

These binary models provide the ability to determine 
the Keplerian parameters of the orbit and various post- 
Keplerian parameters (listed in Table 7). Prom the fitted 
values, TEMP02 can provide various derived quantities in- 
cluding the mass function 

(me sin z)^ _ 477^ (op sin i)^ 
■'~{m^ + mc)2 ~ G P^ ■ ^ ' 

By assuming a typical pulsar mass of rrip — 1.35Mq, a lower 
limit on the companion mass can be estimated by assuming 
that the orbit is viewed edge-on (i = 90°), a median mass 
{i = 60°) and an upper bound at the 90% confidence level 
from i = 26.0° (see, e.g., Lorimer & Kramer 2005). These 
values are obtained by solving the mass function for the com- 
panion mass, rric, using a Newton- Raphson method. If sini 
and Trie are known (e.g. from Shapiro delay measurements) 
then the mass function gives the pulsar mass. 
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Table 6. Astrometric parameters included in tempo2. 



Parameter 


Description 


Symbol 


RA 


Right ascension of pulsar (hr min sec) 


a 


DEC 


Declination of pulsar (deg min sec) 


5 


ELONG 


Ecliptic longitude of pulsar (deg) 


A 


ELAT 


Ecliptic latitude of pulsar (deg) 


P 


PMRA 


Proper motion in right ascension (masyr^-'^) 


Ma 


PMDEC 


Proper motion in declination (masyr^^) 


^^5 


PMELONG 


Proper motion in ecliptic longitude (masyr~^) 


MA 


PMELAT 


Proper motion in ecliptic latitude (masyr"'^) 


M/3 


PX 


Parallax (mas) 


n 


PMRV 


Radial proper motion 


Mil 


POSEPOCH 


Position epoch (MJD) 


ipos 



Table 7. Binary orbital parameters included in tempo2. 



Parameter Description 



Symbol 



PB Binary period of pulsars (d) 

ECC Eccentricity of orbit 

Al Projected semi-major axis of orbit (It-s) 

TO Epoch of pcriastron (MJD) 

OM Longitude of periastron (deg) 

TASC Epoch of ascending node (MJD) 

EPSl esino; 

EPS2 ecoso; 

KOM longitude of the ascensing node 

KIN inclination angle 
SHAPMAX - ln(l - sin i) 



Pb 

e 

X 

To 
ui 

^asc 
K 



OMDOT Pcriastron advance (dcg/yr) 

PBDOT 1st time derivative of binary period 

ECCDOT Rate of change of eccentricity (s^^) 

AIDOT Rate of change of semimajor axis (It-s s~^) 

GAMMA Post-Keplerian 'gamma' term (s) 

XPBDOT Rate of change of orbital period minus GR prediction 

EPSIDOT Rate of change of EPSl 

EPS2DOT Rate of change of EPS2 

MTOT Total system mass (M©) 

M2 Companion mass (Mq) 

DTHETA Relativistic deformation of the orbit 

XOMDOT Rate of periastron advance minus GR prediction (deg yr~^) 

SINI Sine of inclination angle 

DR Relativistic deformation of the orbit 

AO First aberration parameter 

BO Second aberration parameter 

BP Tensor multi-scalar parameter 

BPP Tensor multi-scalar parameter 

AFAC Aberration geometric factor 



UI 

e 

X 

7 

■n 

k 

M 

m2 



s 

dr 

A 
B 

P' 



BPJEP_k Epoch of a step-jump in the binary parameters 

BPJPH_k Size of phase jump 

BTJPB_k Size of jump in orbital period 

BTJAl_k Size of jump in projected semi-major axis 

BTJECC_k Size of jump in orbital eccentricity 

BTJOMJc Size of jump in longitude of periastron 



5.4 Determining pulsar parameters 

Timing models contain the pulsar's rotational, positional 
and binary parameters at a specific epoch. It is often useful 
to determine such parameters at a different epoch. Tempo2 
provides tools to calculate the parameters at any given 



epoch. The position (a, 5) in equatorial coordinates are up- 
dated from the position epoch using proper motion deter- 
minations (/Xa cos (5, Us) 

a = a + ^oL^t a/ cos 5 (16) 
5' = 5 + ^J^sMA (17) 
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where At a represents the difference between the requested 
epoch and the current model position epoch. 

The pulsar's spin frequency and its first derivative are 
updated from the first and subsequent frequency derivatives 



!>' = z> + uAtp + . . 



(18) 
(19) 



where Atp represents the change between the requested 
epoch and the current epoch for the given frequency de- 
terminations. 

Binary parameters are more problematic. The epoch of 
periastron To (or the epoch of the ascending node, Tasc) is 
updated to the closest periastron to the requested epoch by 
calculating the nearest integer n to 



Ah 1 ■ 



Atb 
Pb 



(20) 



which represents the number of orbits since the timing model 
epoch of periastron To, i.e., in interval Atb- The epoch of 
periastron is subsequently updated 



T^ = To + {AtbY 



(21) 



where (Atb)' is obtained by solving Equation 20 with an = 
n. The orbital period (Pb), longitude of periastron (a;), ec- 
centricity (e) and projected semi major axis (oi) are updated 
as 



Pi = Pb + PbiAtbY 

U>' = U) + L0{Atby 

e =e + e{Atb)' 
a'l =ai-\- di(Atb)' 



5.5 Publishing timing ephemerides 



(22) 
(23) 
(24) 
(25) 



TempoI provided only limited output formats for the tim- 
ing model parameters. The publish plug-in for tempo2 
provides the user with the ability to produce output in a 
table format (as in Table 8 where, in this case, the un- 
certainty on each parameter for PSR J0437— 4715 has been 
multiplied by a factor of two and the error corresponds to 
the uncertainty on the last quoted digit). 

If many parameters are included in the fit, then it is pos- 
sible that the resulting uncertainties are covariant. This can 
be checked using the matrix plug-in that displays the corre- 
lation matrix for the fit (a typical output is shown in Table 9 
for PSR J0437— 4715). It is common for low-eccentricity bi- 
nary pulsars that the epoch of periastron and the longitude 
of periastron are near-degenerate. It is common in publica- 
tions to quote highly covariant parameters with more pre- 
cision than suggested by the formal errors (see for instance 
Ryba & Taylor 1991), however, the use of the ELLl binary 
model is preferred. 

Published timing parameters are used to make predic- 
tions for on-line observations of pulsars, for comparing re- 
sults between different observing systems and for searching 
for variations in the parameters with time. It is therefore es- 
sential that full details of how the timing model was created 
are published alongside the parameters. For any published 



timing parameters it is necessary to 1) indicate whether the 
uncertainties represent 1 or 2a formal errors on the fitted 
parameters and whether a weighted or non-weighted fitting 
procedure was used, 2) specify the Solar System ephemeris 
used, 3) indicate the tempo2 version number and describe 
which of the default tempo2 algorithms have not been used 
and which non-default algorithms included, 4) provide full 
information of any pre- whitening carried out on the dataset, 
5) define the coordinate system used and 6) provide details 
of the clock correction process. The fitting process should be 
iterated until the pre- and post-fit values are identical indi- 
cating that the fit has converged. If a weighted fit is carried 
out then it is necessary, in order to obtain accurate errors on 
the fitted parameters, to ensure that the reduced x'^-value 
of the fit is close to unity. It is also common practice and 
desirable to convert the measured parameters so that the 
period, position and binary epochs refer to the centre of the 
data span. 



6 ANALYSING THE TIMING RESIDUALS 

Any systematic feature in the post-fit timing residuals indi- 
cates that some effect is not being described by the timing 
model. Analysing such features potentially provides impor- 
tant information about various perturbations, including the 
presence of planetary companions, variations in the interstel- 
lar medium, precession of the neutron star or irregularities 
in the pulsar's rotation and spin-down. The timing residuals 
can also provide an indication that the TOAs are affected 
by calibration problems or other instrumental effects. 

The pre- and post-fit timing residuals can be plotted 
in rmmerous ways. For instance, the PLK plug-in allows the 
user to plot the residuals versus parameters such as day, 
observing frequency, binary phase, observation length or the 
parallactic angle. Other plug-ins allow the user to obtain a 
list of the timing residuals, barycentric arrival times, clock 
corrections etc. in a specified tabular format. 

As the true uncertainty on any fitted parameter is the 
combination of the random and systematic errors, it is nec- 
essary to attempt to quantify the effects of the systematic 
errors present in the data. Ryba & Taylor (1991) discuss two 
methods to do this which are implemented in the errors 
plug- in 

(i) it is possible to plot a histogram of the normalized 
postflt residuals (the value of the residual divided by its es- 
timated uncertainty). If this histogram follows a Gaussian 
distribution then it is likely that the data are not signifi- 
cantly affected by systematic effects (Figure 15a). 

(ii) compute averages of consecutive sets 1, 2, 4, 8 and 16 
normalised residuals and plot the standard deviations within 
each group. Random Gaussian measurement errors produce 
deviations oc 1/ VlV where iV is the number of residuals av- 
eraged (Figure 15b). 

Tempo2 plug-ins provide numerous tools for analysing 
the post-fit timing residuals or for outputting the residuals in 
formats that are suitable for other data-analysis packages. 
We have noted that large numbers of packages have been 
created to analyse timing residuals from the output files of 
the tempo 1 software and therefore provide facilities within 
TEMP02 to produce output files with the same format as 
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Table 8. Example pajrameters in a I^TJjjX table format for PSR J0437— 4715 obtained using the publish plug-in. 
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Fit and data-set 



Pulsar name J0437-4715 

MJD range 53041.3—53767.3 

Number of TO As 603 

Rms timing residual (fis) 1.6 

Weighted fit N 



Measured Quantities 

Right ascension 04:37:15.78858(13) 

Declination -47:15:08.4685(15) 

Pulse frequency (s"!) 173.68794630602(7) 

First derivative of pulse frequency (s~^) —1.7292(4) X 10~^^ 

Dispersion measure (cm-3pc) 2.64123(17) 

Proper motion in right ascension (masyr-i) 120.9(3) 

Proper motion in declination (masyr~^) —71.0(3) 

Parallax (mas) 8(3) 

Orbital period (d) 5.7410464584(16) 

Epoch of periastron (MJD) 51194.620(6) 

Projected semi-major axis of orbit (It-s) 3.36670624(19) 

Longitude of periastron (deg) 0.9(4) 

Orbital eccentricity 0.00001899(11) 



Set Quantities 

Epoch of frequency determination (MJD) 51194 

Epoch of position determination (MJD) 51194 

Epoch of dispersion measure determination (MJD) 51194 

Sine of inclination angle 0.6788 

First derivative of orbital period 3.64xl0~ 

Periastron advance (dc.'K/yr) 0.016 

( 'ouil)ain()ii mass (_V/ . ) {).2'M) 



Derived Quantities 

log j^Q (Characteristic age, yr) 9.20 

log j^g (Surface magnetic field strength, G) 8.76 



Assumptions 

Clock correction procedure TT(TAI) 

Solar system ephemeris model DE405 

Binary model DD 

Model version number 5.00 



Note: Figures in parentheses are twice the nominal la tempo2 uncertainties in the least-significant digits quoted. 



TEMPO 1. However, it is relatively straightforward to develop 
new plug-ins or to convert old software for use by tempo2. 
For instance, a template plug-in is available which users can 
modify for their own specific uses. Many such plug-ins to 
analyse the post-fit timing residuals are currently being de- 
signed or tested. Examples include the implementation of a 
multi-resolution CLEAN deconvolution algorithm and peri- 
odicity searches (M. Rissi; private communication). 



7 PREDICTIVE MODE 

The folding of pulsar signals proceeds on the basis of the 
predicted time evolution of the phase of the pulse train inci- 
dent upon the observatory. The timing model specifies this 
evolution, but is too computationally intensive for reaJ-time 
applications. Like tempoI, tempo2 is able to produce a 



polynomial approximation of the phase, </>(i), and pulse fre- 
quency, z/, over specified time intervals. The number of co- 
efficients and the time span fitted can be set by the user; 
TEMP02 provides a warning message if the rms deviation 
between the model and the data is large and more coeffi- 
cients or a shorter span are necessary. 



7.1 Tempol-compatible polynomials 

To case the transition from tempo 1 to tempo2, a facility 
is provided to produce predictive polynomials in the same 
format as those made by tempo 1. These consist of a series 
of sets of coefficients of polynomials that approximate the 
evolution of pulse phase incident upon the specified obser- 
vatory. Owing to interplanetary and interstellar dispersion, 
the polynomial is specific to a given observing frequency. 
For observations at radio wavelengths, there is usually a sig- 
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Table 9. The correlation matrix for the fitted parameters for PSR J0437— 4715 data obtained using the matrix plug-in for TEMP02. 
The global correlation (gcor) parameter is a measure for the strongest correlation between the fitted variable and a linear combination 
of all other variables, dp = — logj^Q (1 — gcor^)^/^ provides an estimate of the number of "insignificant" digits that should be quoted in a 
timing solution; see text. 



FO 



TO 



Al 



QM 



ECC 



FO 


+1 


00000000 


















TO 


-0 


13474657 


+1 


00000000 














Al 


+0 


11393821 


+0 


00841676 


+1 


00000000 










OM 


-0 


13474630 


+0 


99999998 


+0 


00841640 


+1 


00000000 






ECC 


+0 


00709055 


+0 


18355036 


+0 


10675629 


+0 


18355054 


+1 


00000000 


gcor 


+0 


17835063 


+0 


99999998 


+0 


15566451 


+0 


99999998 


+0 


21251747 


dp 




0.1 




7.8 




0.1 




7.8 




0.1 



(a) 



-2 2 

Normalised residual 



(b) 



^1 2 5 10 20 50 

Number averaged 

Figure 15. a) Normalised residuals for PSR J1909-3714 plotted 

as a histogram with a Gaussian curve overlaid, b) Rms timing 
residuals after averaging different numbers of residuals together 
with a line indicating the expected \/N dependence. These plots 
were created using the errors plug-in. 



nificant variation of instantaneous pulse phase across the 
observing band. If high precision is required, the form of 
this variation is dependent on many parts of the full timing 
model, thereby defeating any gains in simplicity offered by 
substituting a polynomial form for the time evolution. For 
this reason, tempo2 can produce new time- and frequency- 
dependent predictive polynomials (Section 7.2) which are 
recommended for precision applications. 

A simplified approximation to the frequency depen- 
dence may be used in less critical applications. For isolated 
pulsars, the difference in time of emission between radiation 



received simultaneously at an observing frequency, /, and 
the frequency for which the polynomial applies, /o, is: 

An=D (/-^ - /o"') , (26) 

where D is the dispersion constant (Section 3.5). Neglect- 
ing any variation in pulse frequency over the interval Ad 
(safe for isolated pulsars), the resultant phase difference is 
simply i/Ad, where v is the pulse frequency at the epoch of 
consideration. The pulse frequency in the reference frame of 
the observatory is simply obtained from the polynomial, and 
the observing frequency is also generally known by its value 
in the observatory frame. It is therefore convenient to per- 
form the calculation in this frame, after transforming the 
dispersion constant {D, with dimensions of time~^) using 
the "Doppler shift"® value (/3) specified in the polynomial 
data file: 



D 



t, f) ~ 0(t, /o) - u{t, /o) - fo^) — 



(27) 



To the best of our knowledge, the application of f3 to the dis- 
persion constant is often neglected in existing dedispersion 
software, incurring an error of up to ~ Ad • lO^'', typically 
of the order of several microseconds. It should also be noted 
that unlike tempo2, the value of the dispersion measure 
provided in tempo 1 polynomial files does not include the 
interplanetary contribution, typically in the range 100 ns - 
100 /us. 

For isolated pulsars, the limiting factor in the accuracy 
of this approach is the time evolution of the transformation 
between the observatory and barycentric frames (1 + 0). 
Only a single value is provided for each polynomial, which 
typically spans several hours. The rotation of the Earth ac- 
celerates the observatory sufficiently to change 13 at rates of 
up to ~ 10~^° s~^, resulting in errors of the order of tens of 
nanoseconds. Figure 16 illustrates the effect. 

For binary pulsars, more severe effects arc encountered, 
owing to the different orbital phases of emission of radiation 
received simultaneously at different frequencies. The most 
striking result of this is that, for observing frequencies well 
away from the frequency assumed in the construction of the 
polynomial, pulse profile peaks lead or lag the prediction as 
a function of orbital phase. While eye-catching, to the extent 



In TEMP02 the rate correction also includes gravitational red- 
shift and time dilation. 
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Figure 16. Contour map of the difference between the pulse 
phase as predicted by the full tempo2 timing model, versus a 
single-frequency polynomial approximation in conjunction with 

Eq. 27, for an hour-long observation of the binary pulsar PSR, 
J1906-I-0746 with /o = 650 iVIHz. Contours are spaced by 5 ns: 
dashed contours are negative, positive contours are solid and 
the zero contour is dotted. Two main effects are noticeable. The 
largest effect is due to the evolution in the observatory Doppler 
shift, manifested as a monotonic function of time and frequency 
(delay oc t(/~^ — f^'^), where t is the time relative to the obser- 
vation midpoint) . The second, smaller effect is due to second- and 
higher-order terms in the orbital motion of the pulsar (delay ap- 
proximately oc (/~^ — /o~^)^)- At the midpoint in time where the 
first effect vanishes, the second effect appears as a form roughly 
quadratic in frequency offset. Elsewhere it contributes to the over- 
all assymmetry of the difference function. 



Figure 17. Contour map of the fractional difference in apparent 
pulse frequency at a given observing frequency, versus the central 
frequency (650 MIfz), for an hour-long observation of the binary 

pulsar PSR ,11906+0746. Contours are spaced at 10~*, corre- 
sponding to 10 ns of smearing per second of integration. Dashed 
contours are negative. 



persed), Eq. 27 corrects only to first order in Ao. Since the 
topocentric frequency, ^{t, /) itself varies with dispersion de- 
lay, the second order term (expressed as a time delay) is 
given by h.\,a/2c. This effect is apparent in Figure 16, es- 
pecially around t = Q where the differential Doppler effect 
vanishes. The magnitudes of the effect for sample pulsars 
and observing configurations are provided in Table 10. 



that the error in the predicted pulse phase is constant over 
the course of a given integration, this effect should not man- 
ifest directly in the measured pulse arrival times. Likewise, 
if pulse profile data arc combined from several frequencies 
(or an observing band is coherently dedispersed) , as long as 
"dedispersion" proceeds as prescribed in Eq. 27, the appar- 
ent orbital-phase-dependent effect is compensated to first 
order by the variation in the topocentric pulse frequency 
used to convert a dispersion time delay to a phase delay. 
However, in addition to the secular drift of /3, two remain- 
ing effects dictate the use of time- and frequency-dependent 
two-dimensional polynomials for precision timing of binary 
pulsars. 

Firstly, unless taJien into accouut, the variation in ap- 
parent pulse frequency as a function of observing frequency 
can be large enough to cause appreciable smearing of the 
pulse profile if not taken into account. With scintillation, the 
smearing may be significantly asymmetric, resulting in er- 
rors in measured pulse times of arrival. The fractional change 
in pulse frequency is given to first order by Ado/c where a 
is the line-of-sight orbital acceleration. See Figure 17 for an 
example. Smearing widths for sample pulsars and observing 
configurations are listed in Table 10. 

Secondly, if profiles from several frequencies are to be 
combined (or an observing band is to be coherently dedis- 



7.2 Time- and Frequency-dependent predictions 

To overcome the limitations of the tempo 1-type predictive 
polynomial described in the previous section, tempo2 is able 
to compute a two-dimensional polynomial approximation to 
the timing model, with time and observing frequency as its 
arguments. This mode is recommended for precision appli- 
cations in future processing systems. 

The polynomial is expressed in terms of a two- 
dimensional adaptation of the conventional Chebyshev basis 
functions: 

Pij(x,y) = cos(i cos"'' a;) cos(j cos""*" J/). (28) 

Analogous to the one- dimensional discrete form of a 

set of such basis functions is orthogonal if computed on a 

grid of M X N coordinates {xi,yj) satisfying 

cos(A// cos"^ ,T,;) cos(Af cos~^ J/i) = 0. (29) 

The coefficients of the polynomial approximation to a func- 
tion g{x, y) are therefore easily computed via the inner prod- 
uct: 

^ M JV 

Cki = j^^^Pki{xi,yj)g{xi,yj). (30) 
i=i j=i 

This yields an approximation function 
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Table 10. Smearing and phase errors with TEMPOl-style predictive polynomials for sample pulsars 



Name 


P 


DM 


a sin i 


-^orb 




Ai/P 




Aa/P 




(ms) 


(cm pc^'"*) 


(It-s) 


(d) 


(ns) 


(10-6) 


(ns) 


(10-6) 


J1906+0746 


144.1 


217.8 


1.42 


0.17 


10100 


70 


18.8 


0.13 


J0737-3039A 


22.7 


48.9 


1.42 


0.10 


6000 


260 


2.5 


0.11 


B1744-24A 


11.6 


242.2 


0.12 


0.08 


4600 


400 


9.4 


0.8 


B1913+16 


59.0 


168.8 


2.34 


0.32 


3400 


58 


4.9 


0.08 


J1756-2251 


28.5 


121.2 


2.76 


0.32 


2900 


100 


3.0 


0.11 


J1802-2124 


12.6 


149.6 


3.72 


0.70 


1000 


81 


1.3 


0.10 


J1435-6100 


9.3 


113.7 


6.18 


1.35 


345 


37 


0.3 


0.04 


J02 18+4232 


2.3 


61.3 


1.98 


2.03 


26 


11 


14 X 10-3 


6 X 10-3 


B1957+20 


1.6 


29.1 


0.09 


0.38 


16 


10 


4 X 10-3 


2 X 10-3 


J1909-3744 


2.9 


10.4 


1.90 


1.53 


7.6 


2.6 


6 X 10-"' 


2 X 10-"* 


J2145-0750 


16.1 


9.0 


10.16 


6.84 


1.8 


0.11 


14 X 10-5 


8 X 10-6 


B1855+09 


5.4 


13.3 


9.23 


12.33 


0.7 


0.14 


8 X lO-"' 


15 X 10-6 


.IOl:-!7-iT15 


5.8 


2.() 


:-!.:-! 7 


5.71 


0.2 


0.01 


1 X 10 " 


10 



" Maximum smearing due to pulse frequency offset at / = 600 MHz, vs /o = 650 MHz, over 100-second 

integration (see text). Note, values are approximate only: orbital eccentricity neglected. 

^ Maximum timing error due to neglected order-A|j term, for / = 600 MHz, fo = 650 MHz (see text) 



M N 

9{x,y) ^^^CkiPkiix,y), (31) 

fc=i 1=1 

which is exact for x,y coordinates satisfying Eq. 29. For 
computation, this can be rewritten as a one-dimensional 
Chebyshev polynomial with another Chebyshev polynomial 
defining its coefficients: 

N N 

ff(a;,2/)~^TK2/)^Tfc(a;)cfc,, (32) 
1=1 1=1 

where the one-dimensional basis functions, Tk{x) = 
cos(fc cos-^ x), can be computed efficiently using Clenshaw's 
recurrence relation (Press et al. 1992). 

After mapping the requested intervals in time and ob- 
serving frequency to x and y in the interval [—1, 1], tempo2 
computes the coefficients Cki approximating the function 
4>{t,f) — kf^^, where fc is a constant computed to remove 
the bulk of the frequency dependence due to interstellar dis- 
persion. Subtraction of this term significantly reduces the 
number of coefficients needed for an accurate approxima- 
tion. The M X N coefficients and the constant k are written 
to file to allow the later construction of the approximation 
of 4>{t, /) without reference to the full timing model. 

Using this approach, the timing model, including 
its frequency dependence, can be approximated to sub- 
nanosecond accuracy using a modest number of coefficients. 
As an example, the phase evolution of PSR J1906-I-0746 
over the time and frequency interval depicted in Figure 16 
requires nine coefficients in the time axis and six in fre- 
quency, yielding a difference function with an rms variation 
of 450 ps. For an observation of an isolated pulsar, over the 
same time interval and frequency band, only 5x3 coefficients 
are needed to model the time- and frequency-dependent vari- 
ations in phase, which apart from the basic pulsar frequency, 
are dominated by the acceleration of the observatory due to 
the rotation of the earth. 

The software library packaged with tempo2 provides 
routines for evaluating the predicted pulse phase and pulse 



frequency as a function of time and observing frequency. For 
convenience in applications involving folding of pulsar time- 
series data, the software can also produce a piecewise linear 
approximation which minimises the mean error and keeps 
the rms within specified bounds. A plugin called polytest 
is also provided to assess the accuracy of a computed poly- 
nomial (either 1- or 2-dimensional) , reporting on the rms 
and extrema of the residuals and producing plots such as 
those in Figures 16 and 17. 



8 SUMMARY AND CONCLUSION 

In summary we list below the improved features of TEMP02 
compared to tempoI. For coordinate systems, propagation 
delays and ephemerides, tempo2 

• is compliant with lAU 2000 resolutions and implements 

updated precession and nutation models, polar motion, the 
ICRS coordinate system and TCB (SI units) instead of 
TDB; 

• corrects the observing frequency for relativistic time di- 
lation; 

• uses an improved tabulation of the Solar System Ein- 
stein delay; 

• includes atmospheric propagation delays; 

• includes the Shapiro delay due to the planets; and 

• includes the second-order Shapiro delay due to the Sun. 

For the fitting routines, TEMP02 

• has the ability to simultaneously fit to the timing resid- 
uals of multiple pulsars; 

• implements frequency-dependent fitting (not only DM 
delays) ; 

• has the ability to fit for DM at each epoch using simul- 
taneous observations; 

• allows fits for an arbitrary number of pulse frequency 
derivatives, dispersion measure derivatives and glitches; 
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• simplifies and provides flexible methods for placing ar- 
bitrary offsets between TOAs obtained at different frequen- 
cies or observatories; 

• can whiten data using harmonically related sinusoids; 

• provides a brute-force method for obtaining timing so- 
lutions; 

• fully includes the effects of secular motion of the pulsar; 
and 

• includes the orbital parallax terms in binary models. 

The predictive mode in tempo2 provides both time and 
frequency-dependent predictions and therefore deals cor- 
rectly with 

• Earth rotation and 

• binary motion. 

Other miscellaneous improvements include 

• generalised input formats; 

• generalised output formats; 

• the ability to simulate pulse arrival times; 

• the implementation of numerous graphical interfaces; 

• calculating the post-fit residuals (instead of predicting 
them); 

• the ability to update binary parameters to a given epoch 
and 

• generalised clock correction routines. 

We have presented a now software package, TEMP02, 
that supersedes the existing tempo package. Tempo2 has 
been analysed in detail and wc believe that all the correc- 
tions to the measured TOAs that have been implemented are 
accurate to better than 1 ns. The software has been designed 
so that it is easy to modify current routines and to add new 
functions to describe phenomena which affect pulsar timing 
residuals. For instance, a user can easily implement a new 
binary model or create a personal graphical interface. 

The TEMP02 framework is such that it is relatively easy 
to fit global parameters across multiple data sets. This will 
be the topic of a forthcoming paper and will be used in 
the hunt for gravitational waves, refining the Solar System 
ephemeris and establishing a pulsar-based timescale. 
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APPENDIX A: LIST OF PLUG-IN PACKAGES 
AVAILABLE FOR TEMP02 

• Basic, plots a P — P diagram and indicates the po- 
sition of the pulsar being analysed. Options are available 
to also display a sky-projection that indicates the pulsar's 
position and derived parameters such as the pulsar's char- 
acteristic age and surface dipole magnetic field strength are 
determined. 

• Compare, accepts two input parameter files and pro- 
vides routines to compare differences between the residu- 
als obtained using the two models. For instance, this inter- 
face can be used to compare the effects of different clock or 
ephemcris files or different binary models. 

• CompareRes, accepts two input arrival-time files and 
provides routines to compare differences between the resid- 
uals. 

• Delays, shows the tempo2 calculated delays being 
added to the measured arrival times. For instance, clock 
corrections, ephemeris delays and dispersion delays are in- 
cluded. 

• Errors, used to study systematic and random errors 
in the timing residuals. 

• Fake, allows the user to create simulated arrival times 
that fit a given timing model. The addition of red and white- 
noise is possible. 

• General, a user-specified output format for the pre- 
and post-fit parameters. 

• General2, a user-specified output format for display- 
ing the site and barycentric arrival times, the timing resid- 
uals and various clock and propagation corrections. 

• Gorilla, finds timing solutions using a brute-force 
method. 

• List, provides a listing of the arrival times, residuals, 
clock corrections and propagation delays. 

• Matrix, displays the correlation matrix of the fitted 
parameters. 

• Plk, an interface that plots the timing residuals versus 
parameters such as day, binary phase, length of observation 
etc. Various functions are available including the ability to 
highlight selected points, view pulse profiles and to delete 
observations and refit the data. 

• POLYTEST, provides diagnostics of the approximation 
error of a predictive polynomial, including minimum, max- 
imum and rms and frequency- and time-dependent contour 
maps. 

• Publish, produces publication-quality tables of 
the parameters. 



• Spectral, provides basic spectral analysis tools for the 
timing residuals - including periodograms, auto-correlation 
functions and CLEAN deconvolution. 

• Splk, allows the plotting of multiple pulsar timing 
residuals siirmltancously. 

• Stats, output mode that provides basic information 
about the pulsar and its fit. This plug-in gives the rms resid- 
uals for each observing frequency available. 

• Stridefit, allows fitting of subsets of adjacent obser- 
vations. The resulting fit parameters can be stored and sub- 
sequently plotted versus time. 

• Transform, transforms a timing model created using 
TEMPO 1 to one using the TCB timescale and hence suitable 
for use with tempo2. 



